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Abstract 

In  order  to  apply  a  numerical  analysis  to  design  an  actual  scale  separator  in  polymer  electrolyte  fuel  cell  (PEFC),  it  is  needed  to  enlarge 
the  calculation  size.  In  this  study,  mass  transfer  and  flow  in  gas  diffusion  layer  were  calculated,  and  gas  diffusion  layer  (GDL)  mass  transfer 
approximate  model  based  on  the  theoretical  model  was  developed.  Next,  with  this  model,  PEFC  reaction  and  thermal  flow  analysis  model, 
which  enabled  us  to  calculate  the  case  of  an  actual  scale  cell,  was  made.  Furthermore,  the  effects  of  separator  channel  depth  on  output 
performance  and  current  density  distribution  were  examined  with  this  numerical  analysis.  As  a  result,  in  the  case  of  shallow  channels,  the 
oxygen  transfer  rate  to  electrode  increased  because  of  gas  flow  in  GDL.  However,  current  density  distribution  and  pressure  drop  increased, 
too.  This  calculation  model  can  help  us  to  design  the  optimal  separator  shape  from  the  comprehensive  viewpoint. 

©  2005  Elsevier  B  .V.  All  rights  reserved. 
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1.  Introduction 

Recently,  humankind  has  serious  problems  of  the  environ¬ 
ment  and  energies,  for  example,  global  warming,  acid  rain  or 
lack  of  fuel  sources.  In  order  to  contribute  to  solve  these  prob¬ 
lems,  fuel  cell  is  expected  to  be  practical  use  because  it  has 
low  emission  of  the  environmental  pollutant  and  high  con¬ 
version  efficiency  from  chemical  energy  to  electrical  energy. 
Especially,  polymer  electrolyte  fuel  cell  (PEFC)  is  expected 
as  driving  power  of  vehicles  and  stationary  power  supply.  And 
the  performance  of  PEFC  is  greatly  improved  because  of  the 
development  of  new  component  and  optimization  of  the  sys¬ 
tem.  However,  PEFC  is  demanded  to  have  long  life  and  high 
durability  to  be  generalized  more  widely  as  soon  as  possible. 
The  PEFC  power  generation  characteristic  is  affected  by  the 
structure,  the  material  and  operating  conditions.  The  phenom¬ 
ena  of  mass  transfer,  heat  transfer,  catalysis,  electrochemical 
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reaction  and  fluid  dynamics  are  shown  in  the  internal  cell, 
and  it  is  greatly  important  to  understand  the  correlation  of 
such  complex  phenomena  in  detail  to  improve  and  optimize 
the  PEFC  component  and  system.  However,  regardless  of  the 
size  of  area  from  the  interface  of  catalyst  layer  and  the  stack, 
these  phenomena  are  caused.  And  these  phenomena  affect 
each  other  intricately.  Therefore,  it  is  very  difficult  to  mea¬ 
sure  local  condition  accurately  in  experiments,  and  very  few 
researchers  examine  these  things. 

Recently,  a  numerical  analysis  method  has  been  used  to 
examine  them.  Bemardi  and  Verbrugge  [1,2]  and  Springer 
et  al.  [3]  developed  a  one-dimensional  model  to  the  direc¬ 
tion  of  membrane  thickness,  and  examined  concentration 
distribution  and  water  management  in  PEFC.  Fuller  and  New¬ 
man  [4]  developed  a  two-dimensional  model  to  the  direction 
of  membrane  thickness  and  gas  flow  channel.  Nguyen  and 
White  [5]  and  Yi  and  Nguyen  [6]  developed  heat  and  water 
transport  models  (2D)  that  accounted  for  various  operation 
conditions  and  membrane  hydration  conditions.  On  the  other 
hand,  it  is  thought  that  analysis  with  the  computational  fluid 
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CJ\ 

heat  flux  from  gas  phase  to  solid  phase 
(J  m-2  s-1) 

615 

heat  flux  from  back  plate  to  solid  phase 
(J  m-2  s-1) 

616 

latent  heat  value  of  condensation  (J  m-2  s-1) 

Qb 

all  gas  flow  rate  through  GDL  per  unit  volume 
to  next  channel  (s-1) 

Qb(n) 

flow  rate  through  GDL  per  unit  volume  to  next 
channel  of  n  direction  (s-1) 

Qb(n,m) 

inflow  rate  through  GDL  per  unit  volume  from 
next  channel  of  n  direction  (s-1) 

Qb(n,out) 

outflow  rate  through  GDL  per  unit  volume  to 
next  channel  of  n  direction  (s-1) 

O' 

molar  flux  of  species  j  (molm-2  s-1) 

R 

gas  constant  (8.3 14  J  mol-1  K-1) 

^ohm 

resistance  of  proton  transfer  through  elec¬ 
trolyte  membrane  (Q  m2) 

^rea 

all  reaction  rate  (s-1) 

Re 

Reynolds  number  defined  in  Eq.  (8) 

Sc 

Schmidt  number  defined  in  Eq.  (9) 

Sh 

Sherwood  number  defined  in  Eq.  (20) 

t 

time  (s) 

T 

gas  phase  temperature  (K) 

Tn 

gas  temperature  in  next  channel  of  n  direction 
(K) 

r 

solid  phase  temperature  (K) 

jb 

back  plate  temperature  (K) 

U 

average  gas  velocity  in  GDL  of  v  direction 
(ms-1) 

Uj 

overall  heat  transfer  coefficient  between  gas 
and  back  plate  (J  m-2  s-1  K-1) 

overall  heat  transfer  coefficient  between  back 
plate  and  solid  phase  (J  m-2  s-1  K-1) 

V 

flow  velocity  (ms-1) 

V 

operation  voltage  (V) 

wc 

channel  width  (m) 

wL 

land  width  (m) 

X 

distance  in  v  direction  (m) 

y 

distance  in  y  direction  (m) 

Greek  letters 

a 

net  water  transfer  coefficient 

at 

transfer  coefficient 

p 

fitting  parameter  defined  in  Eq.  (11) 

£ 

effective  porosity  of  GDL 

y 

variable  defined  in  Eq.  (7)  (A mmol-1) 

X 

parameter  defined  in  Eq.  (21) 

i± 

viscosity  of  mixture  gas  (Pas1) 

p 

density  of  mixture  gas  (kgm-3) 

Superscripts 

a 

anode 

c 

cathode 

Nomenclature 

bc 

condensation  rate  constant  (s-1) 

Cj 

molar  concentration  of  species  j  (mol  m-3) 

Cj(n) 

molar  concentration  of  species  j  in  next  chan¬ 
nel  of  n  direction  (mol  m-3) 

re 

co2 

oxygen  concentration  at  catalyst  layer 
(molm-3) 

co2 

oxygen  concentration  at  v  =  0  in  Fig.  7 
(molm-3) 

z^ref 

co2 

reference  oxygen  concentration  (mol  m-3) 

specific  heat  at  constant  pressure  (J  kg-1  K-1) 

P 

diffusion  coefficient  of  species  j  (m2  s-1) 

Df 

effective  diffusion  coefficient  of  species  j 
(m2  s-1) 

erf 

error  function 

E 

electromotive  force  (V) 

Eah 

the  value  of  reduction  change  of  water  enthalpy 
to  voltage  (V) 

F 

Faraday’s  constant  (96,485  C  mol-1) 

h 

heat  transfer  coefficient  of  gas  (Jm-2  s-1  K-1) 

#GDL 

length  of  GDL  gas  flow  area  in  Fig.  7  (m) 

A//H2o 

change  of  water  enthalpy  between  vapor  and 
liquid  (J  mol-1) 

i 

current  density  (Am-2) 

io2 

oxygen  exchange  current  density  (Am-2) 

k 

thermal  conductivity  of  solid  phase 
(Jin-1  s_1  K_1) 

kp 

permeability  of  GDL  (m2) 

kSQ  P 

thermal  conductivity  of  separator 

(Jm-1  s_1  K_1) 

^d,g 

gas  channel  depth  (m) 

Igdl 

GDL  thickness  (m) 

Is 

thickness  of  solid  phase  (m) 

/se  p 

separator  thickness  between  back  plate  and  gas 
phase  (m) 

Mj 

molecular  weight  of  species  j  (kg  mol-1) 

^o2  (reaction)  oxygen  reaction  rate  at  electrode 
(molm-2  s-1) 

/ve 

7V0  2(Re) 

increased  oxygen  mole  flux  to  electrode  by  gas 
flow  (molm-2  s- 1 ) 

/Ve 

yvO  2(Re= 

0)  oxygen  mole  flux  to  electrode  at  Re- 0 
(molm-2  s-1) 

P 

pressure  in  Eqs.  (2),  (3),  (28)  and  (36)  (Pa) 

Pn 

pressure  in  next  channel  of  n  direction  (Pa) 

^H20,sat 

saturated  vapor  pressure  in  stream  (Pa) 

Pe 

Peclet  number  defined  in  Eq.  (23) 

q\ 

heat  flux  from  solid  phase  to  gas  phase 
(J  m-2  s-1) 

612 

heat  flux  from  back  plate  to  gas  phase 
(J  m-2  s-1) 

613 

heat  value  generated  by  reaction  (J  m-2  s-1) 
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channel 

channel 

e 

electrode 

eff 

effective 

k 

anode  or  cathode 

s 

solid  phase 

sep 

separator 

Subscripts 

ave 

average 

h2o 

water 

h2o(1) 

liquid  water 

H20(v) 

vapor  water 

j 

species  j 

n2 

nitrogen 

02 

oxygen 

V 

v  direction 

J 

y  direction 

dynamics  (CFD)  technique  is  important  to  calculate  the  trans¬ 
port  phenomena  in  detail,  and  studies  of  this  kind  have  been 
increasing  recently.  Um  et  al.  [7]  and  Wang  et  al.  [8]  devel¬ 
oped  the  two-dimensional  model  with  CFD,  which  included 
two-phase  flow.  Dutta  et  al.  [9]  made  a  three-dimensional 
computational  model  based  on  a  commercial  software  pack¬ 
age  (Fluent).  Berning  et  al.  [10]  presented  a  non-isothermal 
and  the  three-dimensional  models,  and  calculated  the  distri¬ 
bution  of  current  density  and  concentration  in  straight  chan¬ 
nels.  Mazumder  and  Cole  [11]  examined  the  liquid  water 
transport  with  the  three-dimensional  model.  Li  et  al.  [12] 
analyzed  the  flow  and  concentration  distribution  in  a  small 
cell  with  the  three-dimensional  analysis.  Berning  and  Djilali 
[13]  examined  the  effect  of  porosity  and  thickness  of  gas  dif¬ 
fusion  layer  in  straight  channels  with  the  three-dimensional 
model.  These  PEFC  numerical  analysis  models  contributed 
to  the  optimization  of  component  design  and  operation  con¬ 
ditions,  and  the  examination  of  issues  included  in  the  present 
cell. 

In  PEFC,  anode  and  cathode  gases  usually  flow  through 
each  channel.  And  the  reactant  gas  diffused  to  the  interface 
of  membrane  electrode  assembly  (ME A)  through  gas  dif¬ 
fusion  layer  (GDL).  GDL  is  porous,  and  it  plays  a  role  to 
help  hydrogen  and  oxygen  move  to  an  electrode  catalyst, 
to  support  MEA  and  to  pass  the  electron.  The  ordinary  gas 
flow  situation  is  shown  in  Fig.  1.  However,  in  the  case  of 
a  large-scale  cell,  it  is  thought  that  the  differential  pressure 
between  adjoining  channels  increases,  and  that  the  supplied 
gas  flows  through  GDL  owing  to  the  differential  pressure. 
This  is  shown  in  Fig.  2.  As  a  result,  it  is  supposed  that 
the  cell  endurance  decreases  because  the  temperature  and 
the  humidity  condition  are  not  uniform.  However,  it  can  be 
expected  that  oxygen  transfer  rate  to  electrode  and  output 
density  increase  by  such  gas  flow  in  GDL.  Therefore,  it  is 
important  to  examine  the  influence  of  this  gas  flow  on  cell 


Gas  flow 

along  the  channel 


Diffusion 


Fig.  1.  Ordinary  gas  flow  condition  in  PEFC. 


performance  and  to  optimize  operating  condition  and  cell 
shape. 

Nguyen  [14]  proposed  the  interdigitated  channel  that  had 
the  closed  channel,  and  supplied  gas  was  forced  to  flow 
through  GDL.  Um  and  Wang  [15]  examined  gas  flow  in 
this  interdigitated  channel  with  two  phases  and  the  three- 
dimensional  model.  But,  in  the  case  of  a  large-scale  cell, 
there  was  a  possibility  that  such  gas  flow  through  GDL  was 
occurred  by  large  pressure  drop  in  usual  channel  shape.  Dohle 
et  al.  [16]  and  Oosthuizen  et  al.  [17]  mentioned  this  gas  flow 
through  GDL  with  serpentine  channel  and  examined  gas  flow 
rate  distribution  experimentally  and  numerically.  But  current 
density  distribution  and  cell  performance  were  not  examined. 

In  our  past  researches  [18],  the  effects  of  changing  oper¬ 
ation  temperature,  humidify  temperature  and  hydrogen  and 
oxygen  concentration  in  supply  gas  on  the  i-V  characteristic 
of  a  small  PEFC  were  examined  experimentally.  For  the 
experiment,  we  developed  two  models:  one  was  PEFC 
reaction  model  that  could  show  these  influences  on  PEFC 
reaction  characteristics,  and  the  other  was  PEFC  reaction  and 
flow  analysis  model  that  was  combined  with  the  thermal  flow 
analysis.  With  this  PEFC  reaction  and  flow  analysis  model, 
five  kinds  of  separators  were  evaluated  from  the  viewpoint 
as  follows:  gas  flow  condition,  uniformity  of  current  density 
and  temperature,  reduction  of  pressure  drop  and  ejection  of 
water.  In  [19],  the  simplified  two-dimensional  PEFC  analysis 


Gas  flow 

along  the  channel 


Fig.  2.  Gas  flow  through  gas  diffusion  layer. 
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model  including  flow  and  heat  transfer  of  cooling  water  was 
made.  And  the  influence  of  changing  thickness  of  membrane 
and  GDL  on  cell  performance  was  calculated.  However, 
these  our  past  models  did  not  include  the  effect  of  gas  flow 
through  gas  diffusion  layer.  As  the  gas  flow  rate  and  concen¬ 
tration  and  temperature  are  not  uniform  in  an  actual  cell,  it  is 
also  important  to  examine  the  effect  of  internal  phenomena 
on  cell  performance  in  an  actual- sized  cell  with  an  area  of 
10  cm2.  This  examination  enables  to  identify  the  factors 
which  influence  cell  output  performance  and  durability,  and 
to  optimize  the  PEFC  system.  In  this  study,  especially  gas 
flow  through  GDL  was  aimed,  and  our  former  analysis  model 
was  improved.  The  effect  of  gas  flow  through  GDL  on  the  cell 
performance  and  internal  phenomena  of  an  actual  PELC  sin¬ 
gle  cell  was  examined.  As  mentioned  above,  in  other  studies 
which  examined  concentration  and  liquid  water  distribution 
in  GDL  by  a  numerical  analysis  including  gas  flow  through 
GDL,  the  strict  analysis  model  based  on  CFD  with  an  area 
of  several  square  millimeters  was  mainly  used.  Therefore, 
from  the  viewpoint  of  calculation  resources  and  calculation 
time,  it  is  difficult  to  extend  this  former  model,  which  was 
used  in  other  studies  to  an  actual- sized  cell.  In  this  study, 
the  numerical  analysis  model  was  separated  two  models  in 
consideration  of  practical  calculation  resources  and  practical 
calculation  time.  First,  two-dimensional  mass  transfer  and 
flow  in  cathode  GDL  were  calculated,  and  the  effect  of  gas 
flow  through  GDL  on  oxygen  mass  transfer  rate  to  electrode 
was  examined  under  various  conditions.  This  results  and 
theoretical  mass  transfer  model  were  combined,  and  the 
approximate  model  of  oxygen  mass  transfer  in  cathode  GDL 
was  developed.  Second,  with  this  model,  the  quasi-two- 
dimensional  PEFC  reaction  and  thermal  flow  analysis  model, 
which  enabled  to  calculate  an  actual- sized  cell  was  made. 
Furthermore,  the  effect  of  separator  channel  depth  on  output 
performance  and  current  density  distribution  was  examined 
by  this  numerical  analysis  including  gas  flow  though  GDL. 


2.  Numerical  analysis  model  including  gas  flow 
through  GDL 

In  this  study,  the  PEFC  numerical  analysis  model  was 
developed  step  by  step.  First,  equations  of  motion  and  mass 
balance  in  gas  diffusion  layer  were  calculated  by  the  finite  dif¬ 
ference  method.  And  GDL  mass  transfer  approximate  model 


including  the  effect  of  gas  flow  through  GDL  was  made.  Next, 
with  this  approximate  model,  PEFC  reaction  and  thermal  flow 
analysis  model,  which  could  calculate  an  actual  scale  cell, 
was  developed. 


2.7.  Mass  transfer  and  flow  analysis  in  GDL  and 
development  of  approximate  model 


PEFC  consists  of  ME  A,  two  GDLs  and  two  separators. 
Polymer  membrane  of  MEA  has  proton  conductivity  and 
is  sandwiched  between  two  thin  platinum  electrode  layers. 
MEA  is  sandwiched  between  two  GDLs  and  two  separators. 
In  this  study,  it  was  assumed  that  hydrogen  transfer  rate  in 
anode  GDL  was  faster  than  the  other  mass  transfer  rate  and 
reaction  rate  in  PEFC,  and  only  oxygen  transfer  in  cathode 
GDL  was  examined.  Fig.  3  shows  a  schematic  of  mass  trans¬ 
fer  and  flow  analysis  in  cathode  GDL.  Thickness  of  GDL 
(7gdl)  was  300  pm,  and  width  of  channel  (wq)  and  land  (w\f) 
were  both  1  mm.  The  governing  equations  in  this  simulation 
were  derived  from  the  following  assumptions: 


1 .  The  effective  porosity  and  the  permeability  were  uniform 
in  GDL. 

2.  Water  condensation  was  ignored  in  GDL. 

3.  Density,  viscosity  and  diffusion  coefficient  of  cathode  gas 
were  not  uniform  actually  because  the  composition  was 
changed  locally.  However,  in  this  study,  these  values  were 
regarded  as  constant  and  uniform  for  convenience. 

4.  Electrolyte  membrane  was  humidified  well,  and  ionic  con¬ 
ductivity  was  constant  and  uniform. 

5.  Water  transfer  and  gas  crossover  through  membrane  were 
ignored. 

6.  Gas  flow  condition  in  GDL  was  laminar  flow. 

7.  Cell  voltage  was  uniform. 

8.  The  difference  of  concentration  between  adjoining  chan¬ 
nels  was  very  little,  and  the  concentration  of  channels  was 
set  to  be  equal  to  each  other. 

The  equation  of  continuity  is  shown  by  the  following 

expression: 


fyVx  fyVy 
dx  3  y 


where  p  is  density  of  mixed  gas  and  vx  and  vy  are  velocity  of 
v  and  y  directions. 
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Fig.  3.  Calculation  model  of  cathode  gas  diffusion  layer  in  PEFC. 
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The  equation  of  motion  in  gas  diffusion  layer,  which  was 
porous  media,  was  derived  from  Darcy’s  model. 


kp  dp 

Vx  =  Q 

pi  dx 

(2) 

kp  dp 

Vy=  „ 

r  dy 

(3) 

where  kp  is  the  permeability  of  GDL,  /x  the  viscosity  of  mixed 
gas  and  p  is  the  pressure. 

The  equation  of  mass  balance  of  cathode  gas  is  shown  by 
the  following  expression: 


a  pff  d2Ci  Pffd2cj 

(C  jVy)  +  Df  -r-J  +  Df  3 


3 y 


dx2 


J 


dy- 


(4) 


where  Cj  is  the  concentration  of  chemical  species  j  and  Dy 
is  the  effective  diffusion  coefficient  of  j.  This  equation  was 
derived  to  oxygen,  vapor  and  nitrogen.  The  effective  diffu¬ 
sion  coefficient  was  calculated  with  effective  porosity  by  the 
following  equation: 

Df  =  sDj  (5) 

In  this  study,  overvoltage  of  anode  reaction  was  ignored, 
and  the  relationship  between  cell  voltage  and  overvoltage  is 
shown  by  the  following  equation: 


RT 

at2F 


^ohncd 


where  E  is  electromotive  force,  R  the  gas  constant,  T  the  tem¬ 
perature,  at  the  transfer  coefficient,  F  the  Faraday’s  constant, 
i  the  current  density,  io2  the  exchange  current  density,  Cq^ 
the  reference  oxygen  concentration,  Cqo  the  oxygen  concen¬ 
tration  at  an  interface  of  catalyst  layer  and  Rohm  is  resistance 
of  proton  transfer  through  the  electrolyte  membrane.  R0hm 
was  calculated  by  Nguyen’s  equation  [5].  Reference  oxygen 
concentration  and  exchange  current  density  in  Eq.  (6)  were 
shown  as  the  variable  y  by  the  following  equation: 


*o2 

z^ref 

co2 


In  order  to  calculate  Eqs.  (l)-(6),  boundary  conditions 
were  given  by  the  following  equations: 


(Face  A-B,  C-D  and  E-F) 

dp  _  dCo2  _  dCu20  _  9Gn2 

dy  ’  dy  dy 

(Face  B-C  and  D-E) 

PBC  —  Pin i  PDE  —  Pout 

_  /^channel  _  ^channel 

ch2o  —  ch2o  ’  cN2  —  <-n2 


dy 


r'  _  ^channel 

co2  —  e02 


Table  1 


Calculation  condition  of  GDL  mass  transfer  analysis 


Differential  pressure  between  channels  (Pa) 

0-1000 

Temperature  (°C) 

60 

Transfer  coefficient 

0.3 

Effective  porosity 

0. 1-0.3 

Cathode  gas  density  (kgm-3) 

0.977 

Cathode  gas  viscosity  (Pa  s) 

1.81  x  10"5 

Permeability  of  GDL  (m2) 

2.50x10“ 11 

Electromotive  force  (V) 

1.23 

Membrane  thickness  (fxm) 

30 

Oxygen  diffusion  coefficient  (m2  s-1) 

2.00  x  10“5 

y  of  Eq.  (7)  (A  m  mol- 1 ) 

4.0  x  10“2 

(Face  G-H) 


dp  pi  —i 
dy  kp  4  Fp 

dCu2Q  _ _ i_ 

dy 


Mq2  —  2Mh2o] 


3Cn; 


dCp- 
d  y 


i 


4FDf  ’ 


2^<o 


(Face  A-G  and  F-H) 


dy 


=  0 


dp 


=  0, 


dx 

Cj  is  periodic  boundary  condition  between  A-G  and  F-H 


It  is  necessary  to  note  that  the  analysis  model  of  mass 
transfer  and  flow  in  GDL  was  developed  as  a  local  model  for 
making  the  approximate  model  of  oxygen  mass  transfer  in 
cathode  GDL,  and  that  it  was  not  developed  for  the  examina¬ 
tion  of  the  whole  actual  cell.  In  this  study,  it  was  assumed  that 
periodic  boundary  condition  was  applied  on  the  right  and  left 
sides  of  GDL  (Face  A-G  and  F-H)  in  Fig.  3,  though  it  is  not 
appropriate  in  consideration  of  the  actual  cell  realistically. 

Eqs.  (l)-(6)  were  calculated  with  above  boundary  condi¬ 
tion  by  the  finite  difference  method.  Calculation  condition 
was  shown  in  Table  1.  In  this  study,  the  effect  of  liquid 
water  was  not  directly  included  in  this  analysis  model,  in 
other  words,  the  equations  of  motion  and  mass  balance  of 
liquid  water  were  not  calculated.  So  the  effect  of  diffusion 
inhibition  of  liquid  water  in  GDL  cannot  be  directly  eval¬ 
uated.  However,  in  this  model,  GDL  effective  porosity  was 
used  as  one  parameter  when  the  oxygen  mass  transfer  rate  in 
GDL  was  calculated.  This  value  was  different  from  real  GDL 
porosity,  and  it  was  used  as  apparent  GDL  porosity  includ¬ 
ing  liquid  water.  So  GDL  effective  porosity  was  assumed  as 
0. 1-0.3  to  include  the  effect  of  liquid  water  in  GDL  for  con¬ 
venience,  though  real  GDL  porosity  was  about  0.70.  In  this 
study,  current  density  and  oxygen  concentration  distribution 
were  examined  with  this  analysis  model  changing  differen¬ 
tial  pressure  and  effective  porosity.  And  in  order  to  estimate 
the  effect  of  flow,  dimensionless  numbers  were  defined  by 
the  following  equations: 


Re  = 


iGVhpU 


pi 
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Re  =  3.5 
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Fig.  4.  (a-d)  Oxygen  mole  fraction  distribution  in  GDL  at  0.6  V:  effective  porosity  is  0.3. 


Re  is  the  Reynolds  number,  U  the  average  velocity  of  the  v 
direction  at  x  =  2  mm  and  Sc  is  the  Schmidt  number. 

Figs.  4  and  5  show  oxygen  mole  fraction  distribution  at 
0.6  V  and  0.5  V,  respectively.  In  these  graphs,  in  the  case  of 
(a)  Re  =  0.0,  which  differential  pressure  between  channels 
was  almost  zero  and  gas  did  not  flow,  it  was  found  that  oxy¬ 
gen  moved  to  the  upper  area  of  the  channel  section  more  than 
that  of  the  land  section.  And  comparing  Figs.  4  and  5,  oxygen 
concentration  distribution  became  large  with  low  cell  voltage. 
Next,  changing  the  differential  pressure  between  channels,  it 
was  found  that  oxygen  moved  to  the  upper  area  of  the  land 
section  by  gas  flow,  and  that  the  effect  became  large  accord¬ 
ing  to  Reynolds  number  in  Figs.  4b-d  and  5b-d.  Fig.  6  shows 
current  density  distribution  at  0.6  V.  As  gas  flow  velocity  was 
large,  current  density  was  high  at  the  upper  area  of  the  land 
section.  On  the  other  hand,  current  density  at  the  upper  area 
of  a  downstream  channel  (x  =  2. 0-4.0  mm)  was  affected  by 
convection  flowing  out  to  the  downstream  channel,  and  it  was 
locally  lower  than  that  without  gas  flow.  Moreover,  when  Re 
was  0.0,  1.4,  3.5  and  7.0,  the  average  current  densities  were 
0.862,  0.884,  0.900  and  0.911  A  cm-2,  respectively,  and  it 
was  found  that  the  average  current  density  was  not  propor¬ 
tional  to  the  Reynolds  number.  At  face  A-G  and  F-H,  the 


Fig.  6.  The  effect  of  gas  flow  through  GDL  on  current  density  distribution 
at  0.6  V:  effective  porosity  is  0.3. 

periodic  boundary  condition  was  applied.  However  in  this 
calculation  results,  the  concentration  gradient  at  this  bound¬ 
ary  was  almost  zero,  and  the  mass  flux  through  this  boundary 
by  diffusion  could  be  ignored.  In  other  words,  the  calculation 
results  equaled  to  the  results  with  the  boundary  condition  that 
concentration  gradient  was  zero,  and  this  boundary  condition 
did  not  affect  the  whole  calculation  results. 
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Fig.  5.  (a-d)  Oxygen  mole  fraction  distribution  in  GDL  at  0.5  V:  effective  porosity  is  0.3. 
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Catalyst  layer 


Fig.  7.  Schematic  of  GDL  approximate  model. 


Mass  transfer  with  gas  flow  was  examined  by  the  numer¬ 
ical  analysis  mentioned  above.  However,  it  was  difficult  to 
apply  this  model  to  an  actual  scale  cell  from  the  viewpoint 
of  calculation  time  and  calculation  resource.  So  the  simpli¬ 
fied  mass  transfer  model,  which  expressed  oxygen  transfer 
with  gas  flow  through  GDL,  was  developed.  In  this  study, 
it  was  based  on  a  theoretical  model,  and  given  by  fitting 
with  numerical  simulations.  In  Figs.  4  and  5,  the  section 
from  v  =  0  mm  to  2  mm  was  defined  as  an  upstream  area, 
and  the  section  from  x  =  2mm  to  4  mm  was  defined  as  a 
downstream  area.  The  upstream  area  and  the  downstream 
area  were  decided  by  the  direction  of  gas  flow  through  GDL 
in  this  study,  and  these  were  not  decided  by  the  position 
of  channels.  And  current  density  and  oxygen  concentra¬ 
tion  at  catalyst  layer  were  averaged  in  each  section  with 
2  mm  width.  The  rate  of  oxygen  consumption  by  electrode 
reaction  was  equal  to  the  summation  of  oxygen  flux  to  cat¬ 
alyst  layer  without  gas  flow  and  oxygen  flux  by  gas  flow 
through  GDL.  This  relationship  is  shown  by  the  following 
equation: 


The  following  convectional  diffusion  equation  was  derived 
from  a  mass  balance  equation  of  oxygen  with  above 
assumption: 


r,dC0 2  _  ne ff92c02 

U~dP  ~  °°1^yr 


(12) 


And  boundary  conditions  were  shown  by  the  following  equa¬ 
tion: 

Co2  =  Cg2,  at*  =  0,  y  >  0 
Co2  =  Co2 ’  at*  >  0,  y  =  0 

For  the  purpose  of  getting  exact  solution  of  Eq.  (12)  with 
above  boundary  condition,  mathematical  technique  was 
required.  The  detailed  solving  method  about  this  partial  dif¬ 
ferential  equation  was  described  in  other  text,  for  example 
reference  [20].  In  this  paper,  the  derivation  of  solution  was 
omitted,  and  the  following  solution  was  obtained: 


(13) 


where  erf  is  error  function,  and  it  is  calculated  by  the  follow¬ 
ing  equation: 


(14) 


The  oxygen  mole  flux  at  interface  of  catalyst  layer  was  cal¬ 
culated  by  the  following  equation: 


^02 (reaction)  ~  ~  NO2(Re=0)  +  N02(Re)  (1^) 

These  oxygen  fluxes  were  modeled,  respectively. 

First,  oxygen  flux  without  gas  flow  was  calculated  by  the 
following  equation: 

^channel  _ 

a  jQ  _  r^eff  O2  02  w  q  / 1  1  ^ 

NQ2(Re= 0)  —  D02  X  ^ 

where  f$  is  the  fitting  parameter  and  it  is  function  of  channel 
and  GDL  shape. 

Next,  oxygen  flux  with  gas  flow  was  examined.  Fig.  7 
shows  the  objective  area  to  develop  the  theoretical  model. 
In  this  figure,  gas  flowed  to  the  *  direction  along  the  cat¬ 
alyst  layer,  but  the  y  direction  needs  to  be  paid  attention 
to  because  it  is  different  from  that  in  Fig.  3.  The  govern¬ 
ing  equations  in  this  model  were  derived  from  the  following 
assumptions: 

1 .  Gas  flow  to  the  y  direction  is  ignored. 

2.  Gas  diffusion  to  the  *  direction  is  ignored. 

3.  Gas  flow  velocity  is  uniform  in  the  y  direction. 

4.  Oxygen  concentration  at  *  =  0  is  constant  and  uniform. 

5.  Oxygen  concentration  at  interface  of  catalyst  layer  (y  =  0) 
is  uniform. 


y=0 


(15) 


Eq.  (16)  was  obtained  by  substituting  Eq.  (13)  for  Eq.  (15) 


Neo2 


(Cm 


o2 


co2) 


Peo2U 

7 XX 


(16) 


The  average  oxygen  mole  flux  of  section  Hqdl  in  Fig.  7  was 
obtained  by  integrating  Eq.  (16). 


^02(ave) 


1 


* Hgdl 


Hgdl  Jo 


Neo2dx 


UPo2 


cb2) 

(17) 


This  average  oxygen  flux  expresses  the  effect  of  gas  flow 
through  GDL.  In  this  study,  Hqdl  was  4  mm,  and  this 
increased  oxygen  flux  was  divided  into  an  upstream  area  and 
a  downstream  area.  And  it  was  assumed  that  Cqo  was  equal 

to  the  average  of  Cq2  and  c^nnel. 


N02(Re)  ~ 


1 

2 


ttHqdl 


/^channel 

lco2 


(18) 
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With  Eqs.  (10),  (11)  and  (18),  the  following  equation  was 
obtained: 


i 

4 F 


/^channel  _ 

co2  co2 


/gdl 


T/qdL 


^channel 


(19) 


The  Sherwood  number  was  defined  by  the  following  equation 
as  a  dimensionless  number,  which  expressed  mass  transfer 
rate: 


Sh 


i  /gdl 

4 F  Df2(C$™nel  -  Cq2) 


(20) 


And  the  variable  X  was  defined  by  the  following  equation: 


x=f  Jgd^ 

2  y  7t//gdl 

Eq.  (19)  was  transformed  with  Eqs.  (20)  and  (21) 
Sh=  P  +  Xs/Ye  =  0  +  xJTeJYc 


(21) 


(22) 


where  Pe  is  the  Peclet  number,  and  it  was  calculated  by  the 
following  equation: 


Pe  = 


UIqdl  UIqdlP  M 


Deff 

u02 


fJ  pD 


eff 

o2 


=  Re  -  Sc 


(23) 


Eq.  (22)  is  a  simplified  mass  transfer  model.  In  order  to  deter¬ 
mine  the  fitting  parameter  /3  and  to  confirm  the  validity  of  this 
equation,  the  relationship  between  Sh  and  Pe  (or  Re  and  Sc) 
was  examined  by  the  numerical  analysis  as  Eqs.  (l)-(7).  The 
average  current  density  (/)  and  the  average  of  difference  of 
oxygen  concentration  (CQh2annel  —  Cq2  )  between  the  upstream 
area  and  the  downstream  area  were  calculated  by  changing 
cell  voltage  from  0.9  to  0.05  with  Eqs.  (l)-(7)  analysis.  Sh 
was  obtained  from  the  gradient  of  those  two  quantities  with 
Eq.  (20).  Moreover,  calculations  of  such  kind  were  repeated 
by  changing  differential  pressure  between  channels  and  by 
changing  effective  porosity. 

Fig.  8  shows  the  relationship  between  Sh  and  Sc  atRe  =  0 
(differential  pressure  was  zero).  It  shows  Sh  is  independent 
of  Sc  and  this  value  is  about  0.624.  From  Eq.  (22),  it  was 
confirmed  that  the  value  of  /3  was  0.624.  Fig.  9  shows  the 
relationship  between  the  Sherwood  number  and  the  Reynolds 
number.  Sh  increased  greatly  by  Re  in  case  of  large  Sc ,  and 
the  mass  transfer  rate  of  an  upstream  area  was  affected  more 
strongly  than  that  of  a  downstream  area  by  Re.  Sh  of  the 
downstream  area  was  not  changed  until  Re  came  to  the  aver¬ 
age  1.3  in  each  Sc,  so  the  quantity  that  1.3  subtracted  from 
Re  was  used  to  examine  Eq.  (22).  Fig.  10  shows  the  rela¬ 
tionship  between  Sh  and  Re -Sc  ( -Pc )  of  the  upstream  area  (a 
and  b)  and  the  downstream  area  (c).  In  Fig.  10a,  a  horizontal 
axis  is  Re-Sc,  and  a  vertical  axis  is  (Sh  —  /3)/X.  If  Eq.  (22) 
is  correct,  all  data  are  expressed  by  a  solid  line.  But  these 
data  were  not  located  on  this  line.  Therefore,  this  equation 


Fig.  8.  Relationship  between  Sherwood  number  and  Schmidt  number: 
Reynolds  number  is  zero. 

was  tried  to  prove  it  correct.  In  Fig.  10b,  the  vertical  axis 
shows  (Sh  —  P)l(XSc 0,3 ).  In  this  graph,  all  data  are  shown  by 
a  solid  line.  It  needs  to  correct  for  the  following  reasons. 
Eq.  (22)  was  derived  from  the  assumption  of  perfect  parallel 
flow  along  the  catalyst  layer.  But  there  was  not  such  flow 
condition  in  the  numerical  simulation,  and  especially  at  the 
upstream  area,  oxygen  in  a  gas  channel  was  moved  to  a  cat¬ 
alyst  layer  by  vertical  gas  flow.  Next  in  Fig.  10c,  in  the  case 
of  the  downstream  area,  the  equation  of  Sh  could  be  modeled 
by  (Re  —  1.3)  in  Eq.  (22).  Consequently,  the  simplified  mass 
transfer  model  as  followings  were  obtained  by  the  theoretical 
models  and  the  numerical  simulations: 

Upstream  Sh  =  0.624  +  XRe°-5Sc0^ 

Downstream  Sh  =  0.624  Re  <  1.3 

Sh  =  0.624  +  X(Re  -  1.3)a55c0-5  Re  >  1.3 

(24) 

Similar  examinations  were  carried  out  as  to  the  different  oxy¬ 
gen  concentration  in  gas  channel  and  ionic  conductivity  of 
electrolyte  membrane,  and  it  was  confirmed  that  the  same 
equation  was  obtained  from  these  examinations.  The  oxygen 
concentration  at  an  interface  of  catalyst  layer  could  be  cal- 


Fig.  9.  Relationship  between  Sherwood  number  and  Reynolds  number. 
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Fig.  10.  Relationship  between  Sherwood  number  and  Reynolds  number- Schmidt  number:  (a)  upstream  area;  (b)  upstream  area  with  correction;  (c)  downstream 
area. 


culated  with  Eq.  (24),  and  the  characteristic  of  cell  voltage 
and  current  density  was  calculated  by  substituting  it  for  Eq. 
(6).  But  this  model  was  given  from  the  assumption  that  liq¬ 
uid  water  was  ignored,  it  has  to  be  improved  in  detail.  As 
the  numerical  value  in  Eq.  (24)  was  obtained  with  the  shape 
of  Fig.  3,  these  values  were  needed  to  examine  in  the  case 
of  other  shapes,  which  meant  different  channel  width,  land 
width  and  thickness  of  GDL.  However,  in  this  study,  the  influ¬ 
ence  of  cell  shape  on  these  values  was  not  examined,  and  the 
calculation  of  the  following  section  was  carried  out  with  the 
same  shape. 

2.2.  PEFC  reaction  and  thermal  flow  analysis  model 

The  oxygen  transfer  model  was  developed  in  the  previous 
section.  In  this  section,  the  PEFC  reaction  and  thermal  flow 
analysis  model,  which  could  calculate  an  actual  scale  cell, 
was  developed.  As  the  definition  of  coordinate  was  different 
from  that  of  previous  section,  and  that  is  an  important  point  to 


notice.  Fig.  11  shows  the  PEFC  simulation  model.  As  shown 
in  this  figure,  gas  flow  velocity,  concentration  and  temper¬ 
ature  were  calculated  in  gas  channels  on  the  anode  and  the 
cathode.  And  it  was  assumed  that  the  temperature  distribu¬ 
tions  of  ME  A  and  GDL  were  the  same  as  each  other  and  they 
were  unified,  and  the  temperature  and  current  density  were 
calculated  in  the  unified  part.  (This  unified  part  is  expressed 
as  a  solid  phase.)  The  governing  equations  in  this  simulation 
were  derived  from  the  following  assumptions: 

1 .  The  inlet  gas  flow  rate  in  each  channel  is  uniform. 

2.  The  volume  of  the  condensation  water  is  ignored,  and 
the  water  moves  with  the  gas. 

3.  Reduction  of  the  reaction  area  caused  by  flooding  of 
electrode  is  ignored,  and  it  is  also  ignored  that  water 
condensation  prevents  the  diffusion. 

4.  Fluid  is  incompressible  Newtonian  fluid  and  ideal  gas. 
Flow  condition  is  laminar  flow.  Gas  properties  are  con¬ 
stant. 
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Fig.  11.  Model  for  the  simulation  of  the  PEFC. 


5.  The  temperature  of  back  plate  is  uniform  and  constant. 

6.  Heat  transfer  between  a  separator  and  gas  is  ignored.  But 
heat  transfer  among  gas  phase,  solid  phase  and  back  plate 
is  included. 

7.  Cell  voltage  is  uniform  and  constant. 

8 .  Only  resistance  overvoltage  and  water  transfer  in  a  mem¬ 
brane  include  the  influence  of  temperature. 

9.  In  a  membrane,  ionic  conductivity,  electro-osmosis  coef¬ 
ficient  and  water  effective  diffusion  coefficient  that 
depend  on  the  membrane  humidity  are  determined  by 
water  activity  of  the  anode  side. 

10.  The  gas  crossover  through  a  membrane  is  disregarded. 

1 1 .  The  permeability  of  GDL  is  constant  and  uniform. 

In  this  study,  the  one-dimensional  analysis  as  plug  flow  in 
each  channel  was  available  because  of  the  assumption  that  the 
inlet  gas  flow  rate  distribution  was  uniform.  Though  the  sepa¬ 
rator  shape  was  a  two-dimension  structure  to  the  direction  of 
the  membrane  faces,  the  quasi-two-dimension  analysis  model 
was  made  by  assuming  the  direction  from  the  inlet  to  the  out¬ 
let  to  be  a  positive  v  direction  in  each  channel  and  by  making  a 
channel  meander.  As  a  result,  simplification  of  the  equations 
and  reduction  in  calculation  time  became  possible.  However, 
in  the  case  of  calculation  of  the  solid  phase  temperature  dis¬ 
tribution,  it  was  calculated  by  the  two-dimensional  analysis 
model.  And  in  order  to  calculate  gas  flow  rate  that  flowed 
to  the  next  channel  through  GDL,  pressure  distribution  was 
calculated  by  the  two-dimensional  analysis.  Moreover,  for 
simplification  of  the  calculation,  these  following  two  terms 
in  gas  channel  were  ignored:  the  heat  conduction  term  in  the 
energy  balance  equations;  and  the  diffusion  term  in  the  mass 
balance  equations. 

The  equation  of  continuity  is  shown  by  the  following  equa¬ 
tion: 

dvk  k  k 

-V  =  -*rea  ~  0b  (25) 


where  v  is  the  velocity  of  mixed  gas,  v  the  distance  along  a 
gas  flow  channel,  7?rea  the  all  reaction  rate,  Q b  the  all  gas  flow 
rate  through  GDL  per  unit  volume  to  the  next  channel  and 
the  superscript  k  is  the  anode  side  or  the  cathode  side.  7?rea  is 
calculated  by  the  following  equation: 


(26) 


where  /d,g  is  depth  of  gas  channel,  p  the  density  of  mixed 
gas,  Mj  the  molecular  weight  of  a  chemical  species  j  and  rj 
is  the  reaction  or  condensation  rate  per  unit  area  of  chemical 
species  j.  Q b  is  calculated  by  the  following  equation: 

0b  =  EeL)  (27> 

n 


where  is  gas  flow  rate  through  GDL  to  the  n  direction. 

The  equation  of  motion  is  shown  by  the  following  equa¬ 
tion: 


-v/  +  Ak(4a  +  gS) 


(28) 


where  p  is  the  pressure,  fi  the  gas  viscosity,  wq  the  width 
of  a  gas  channel  and  the  operator  D/D t  is  substantial  time 
derivative  that  is  shown  by  the  following  equation: 


Duk  dvk  k  dvk 

-  =  - T  v  - 

D t  dt  dx 


(29) 


The  viscous  term  in  Eq.  (28)  is  derived  from  the  Hele-Show 
model.  This  term  includes  the  effect  of  viscous  drags  between 
two  pairs  of  facing  walls  in  a  channel. 

The  equation  of  chemical  species  j  is  shown  by  the  fol¬ 
lowing  equation: 


D  t 


rk 

-A  +  Cik(^rea  +  0b)  +  in) 

n 

-Eck0b(«,ouo  (30) 

n 


where  Cj  is  concentration  of  a  chemical  species  j ,  Cj(n) 
the  concentration  of  a  chemical  species  j  at  the  next  chan¬ 
nel  to  the  n  direction,  Qb(n,m)  the  gas  flow  rate  through 
GDL  from  the  n  direction  adjoining  channel  to  this  point 
and  <2bO,out)  is  gas  flow  rate  through  GDL  from  this  point 
to  the  n  direction  adjoining  channel.  The  equations  of 
species  were  derived  to  eight  kinds  of  chemical  species: 


r^a  /^a  r^a  r^a  /^C  /^C  s~<C  tViQt 

cH2’  cN2’  cH20(v)’  cH20(1)’  c02’  cN2’  cH20(v)’  CH20(1)  lIlcU 
are  hydrogen,  oxygen,  nitrogen,  vapor  and  condensed  water 

in  anode  and  cathode  channel. 
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The  equations  of  energy  are  shown  by  the  following  equa¬ 
tions: 


Reaction  and  condensation  rates  of  each  ingredient  are 
shown  by  the  following  equations: 


Drk  q\  +  cM  k  k  k 

(Gas)  isr  =  +  T +  Cb) 

Dl  p  cpllg 

+ ET«k2b(„,in)  -  E^eU.om)  (3D 

n  n 

(Solid)  psCs  —  =  fcsV2rs  +  ^3  +  <?4  +  *?5  +  <?6  (32) 

p  dt  Is 

In  the  energy  equation  of  gas,  Cp  is  the  specific  heat,  T 
the  temperature  and  q\  and  q2  are  heat  fluxes  from  a  solid 
phase  and  a  back  plate,  respectively.  Tn  is  the  gas  tempera¬ 
ture  in  the  next  channel  to  the  n  direction.  In  the  equation 
of  a  solid  phase,  k  is  the  heat  conductivity,  /s  the  thickness 
of  solid  phase,  <73  the  heat  value  per  unit  area  as  a  result 
of  electrochemical  reaction,  <74  and  q$  the  heat  fluxes  from 
gas  and  a  back  plate,  respectively,  q§  the  latent  heat  flux  of 
water  condensation  and  the  superscript  ‘s’  is  a  solid  phase. 
These  heat  flux  and  heating  value  are  shown  by  the  following 
equations: 


q\  =  hk(Ts  -  Tk) 

q\  =  Uj(Tb  -  Tk) 

q\  =  ( Eah  -  V)i 

q\  =  ha(T a  -  Ts)  +  hc(Tc  -  Ts ) 

qs5  =  uf}(Th  -  Ts )  +  Ujic)(Tb  -  Ts ) 

^6  =  —  A//h2o(^h20(1)  +  rH20(l)) 


(33) 
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where  F is  the  Faraday’s  constant,  R  the  gas  constant,  PH2o,sat 
the  saturated  vapor  pressure,  a  the  water  transfer  coefficient 
and  bc  is  the  condensation  rate  constant. 

Gas  flow  rate  through  GDL  is  calculated  by  the  following 
Darcy’s  model: 


k  /k 
*GDL 


k 

L 


(36) 


where  kp  is  the  permeability  of  GDL,  /qdl  the  thickness  of 
GDL  and  u>l  is  the  width  of  the  land  area  that  is  between 
channels. 

Current  density  i  was  calculated  with  Eqs.  (6)  and  (24).  In 
this  model,  current  density  was  calculated  by  the  following 
equation  as  the  function  of  local  concentration  and  tempera¬ 
ture: 


where  h  is  heat  transfer  coefficient  of  anode  or  cathode 
gas,  Eah  the  value  when  the  change  of  water  enthalpy 
was  converted  to  voltage,  V  the  voltage,  i  the  current 
density,  A//h2o  the  change  of  water  enthalpy  between 
vapor  and  liquid,  Uj  the  overall  heat  transfer  coefficient 
between  anode  gas  or  cathode  gas  and  a  back  plate  and 

Uj  is  the  overall  heat  transfer  coefficient  between  anode 
or  cathode  side  back  plate  and  solid  phase.  These  over¬ 
all  heat  transfer  coefficients  are  shown  by  the  following 
equations: 


u\  = 

Ujk) 


1 


1  l  /sep 

J7  +  fcsip 

fcsep 


/seP  +  / 


k 

4g 


(34) 


^/ffiCii2o(v),C^iJs)  (37) 

It  is  thought  that  water  is  moved  under  two  mechanisms, 
electro-osmosis  and  back  diffusion  in  an  electrolyte  mem¬ 
brane.  When  one  proton  moves  from  the  anode  side  to 
the  cathode  side,  the  water  movement  coefficient  a  shows 
the  net  number  of  water  molecules  moving  along  with 
proton.  This  is  from  the  method  by  Nguyen  and  White 

[5]. 

The  local  concentration,  temperature,  flow  velocity  and 
current  density  were  calculated  with  Eqs.  (25)-(37).  Their 
partial  differential  equations  are  discretized  by  the  finite  dif¬ 
ferential  method.  The  boundary  conditions  of  flow  velocity, 
temperature  and  concentration  are  set  as  followings: 

(1)  The  gas  inlet  boundary :  these  variables  are  constant. 

(2)  The  gas  outlet  boundary :  the  gradients  of  these  variables 
are  constant. 


where  /sep  is  the  separator  thickness  between  a  back  plate  and 
a  gas  phase  and  ksep  is  the  heat  conductivity  of  a  separator. 


Current  density  and  water  transfer  coefficient  were  cal¬ 
culated  all  over  the  electrode  area.  Those  variables  were 
calculated  until  becoming  stationary  state.  The  relative  errors 
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of  balance  equation  of  mass,  species  and  energy  became  less 
than  1%  in  all  the  calculations. 


3.  Results  and  discussions 

3.1.  Relationship  between  gas  flow  rate  and  pressure 
drop  with  or  without  GDL 

In  order  to  confirm  the  validity  of  this  numerical  analy¬ 
sis  model  including  gas  flow  through  GDL,  the  cathode  gas 
flow  rate  and  pressure  drop  were  measured  in  a  25  cm2  cell 
with  or  without  GDL  under  the  condition  without  reactions, 
and  calculation  results  were  compared  with  the  experimental 
results.  Fig.  12  shows  the  gas  channel  shape  of  a  separator  in 
flow  experiments  without  reactions.  Fig.  12a  is  a  serpentine 
channel  separator  with  1  channel,  (b)  is  a  serpentine  channel 


separator  with  5  channels,  (c)  is  a  parallel  channel  separator 
with  25  channels  and  (d)  is  a  semi- serpentine  channel  separa¬ 
tor  with  8  channels.  Width  and  depth  of  the  channel  were  both 
1  mm,  width  of  the  land  area  was  also  1  mm  in  each  separa¬ 
tor,  and  thickness  of  GDL  was  300  pan.  The  supplied  gas  was 
not  humidified,  and  cell  and  gas  temperature  were  25  °C.  The 
value  of  permeability  of  GDL  was  treated  as  2.5  x  10  11  m2 
for  the  numerical  analysis. 

Fig.  13  shows  the  relationship  between  gas  flow  rate  and 
pressure  drop  in  the  cathode  side  of  a  25  cm2  cell  by  the  exper¬ 
iment  and  the  calculation  with  each  separator.  The  influence 
of  gas  flow  through  GDL  was  examined  with  four  kinds  of 
separator  for  the  experiment  and  the  calculation.  However, 
as  this  analysis  model  developed  in  this  study  was  derived 
from  assuming  that  each  channel  did  not  branch  off  or  conflu¬ 
ent  between  the  inlet  and  the  outlet,  the  experimental  results 
could  not  be  compared  with  calculation  results  in  Fig.  13c 


Fig.  12.  (a-d)  Gas  channel  shape  of  separator  in  non-reaction  flow  experiments. 
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Fig.  13.  Relationship  between  gas  flow  rate  and  pressure  drop  in  cathode  side  of  25  cm2  cell  by  experiment  and  calculation:  (a)  serpentine  separator  with  one 
channel;  (b)  serpentine  separator  with  five  channel;  (c)  parallel  separator;  (d)  semi-serpentine  separator. 


and  d.  In  Fig.  13a  and  c,  it  was  found  that  the  influence  of 
GDL  on  pressure  drop  was  large  in  the  case  of  serpentine  sep¬ 
arator  with  one  channel  and  semi-serpentine  separator.  And 
pressure  drop  in  the  experiments  without  GDL  was  1.5  times 
as  large  as  that  with  GDL.  On  the  other  hand,  in  Fig.  13b  and 
d,  the  difference  of  pressure  drop  between  the  experiments 
with  GDL  and  that  without  GDL  was  small.  It  was  considered 
that  the  differential  pressure  between  adjoining  channels  of 
serpentine  separator  with  one  channel  and  semi- serpentine 
separator  was  larger  than  that  of  other  separators.  The  gas 
flow  rate  through  GDL  was  increased  by  its  differential  pres¬ 
sure,  and  the  gas  flow  along  the  channel  became  incomplete. 
Consequently,  the  actual  PEFC  is  sure  to  have  GDL,  which 
occurs  non-uniform  flow,  and  such  flow  condition  is  different 
by  separator  shapes. 

Next,  the  calculation  results  were  compared  with  exper¬ 
imental  results  in  Fig.  13a  and  b.  The  pressure  drop  was 
theoretically  in  proportion  to  gas  flow  rate  when  the  gas  flow 
condition  was  laminar  flow.  However,  as  the  gas  flow  rate  per 
channel  was  high  in  the  case  of  one  channel  separator,  the 
gas  flow  condition  became  turbulence  flow  without  GDL. 
And  pressure  drop  of  the  experiment  was  not  in  proportion 


to  the  gas  flow  rate,  as  a  result,  there  was  a  small  difference 
between  the  experiment  results  and  calculation  results  with¬ 
out  GDL  in  one  channel  serpentine  separator  (Fig.  13a).  In 
other  condition,  the  calculation  results  by  this  model  were 
almost  equal  to  the  experimental  results,  and  the  validity  of 
this  model  including  such  gas  flow  through  GDL  was  con¬ 
firmed. 

3.2.  Effect  of  channel  depth  on  output  performance  and 
current  density  distribution 

Fig.  14  shows  the  separator  shape  of  the  anode  side  and 
the  cathode  side  that  is  the  target  in  this  study.  The  electrode 
area  was  150  mm2.  The  width  of  channel  and  land  was  both 
1  mm,  and  the  number  of  channels  was  15.  Fig.  15  shows  the 
developed  view  of  a  single  cell  and  gas  flow  pattern  with  the 
separator  of  Fig.  14.  As  the  same  separators  were  put  together 
like  Fig.  15a,  gas  flow  pattern  are  shown  in  Fig.  15b.  Table  2 
shows  the  calculation  parameter  including  the  operating  con¬ 
dition  and  the  dimensions  of  MEA,  GDL  and  separator.  It  was 
supposed  that  physical  properties  were  constant,  because  the 
effect  of  changing  gas  composition  in  a  cell  on  physical  prop- 
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Fig.  14.  Gas  channel  shape  of  separator  in  225  cm2  cell  calculation. 

erties  was  little.  GDL  effective  porosity  was  0.2  as  reference 
value,  and  the  effect  of  liquid  water  in  GDL  was  included  in 
analysis  model  indirectly.  And  it  was  supposed  that  the  chan¬ 
nel  depth  of  anode  and  cathode  equaled  each  other.  In  order  to 
examine  the  characteristic  of  current  density-cell  voltage  in 
each  condition,  the  calculation  was  carried  out  from  0.9  V  to 
0.05  V  by  the  0.05  V  steps,  and  it  took  8  h  per  one  calculation 
condition  by  Penitum4®  3.2  GHz  PC. 

The  effect  of  the  channel  depth  was  examined  with  the 
numerical  analysis  model  developed  in  this  study.  Fig.  16 
shows  the  current  density-cell  voltage  curve  of  each  chan¬ 
nel  depth.  In  this  graph,  it  was  found  that  the  cell  voltage  of 
a  shallow  channel  was  higher  than  that  of  a  deep  channel  at 
high  current  density  condition.  The  reason  was  that  the  differ¬ 
ential  pressure  between  adjoining  channels  increased  in  the 
case  of  the  shallow  channel,  and  that  the  oxygen  transfer  rate 
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- ► 
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Table  2 

Calculation  condition  of  PEFC  actual  size  cell  analysis 


Pressure  (MPa) 

0.1 

Inlet  gas  and  humidify  temperature  (°C) 

60 

Back  plate  temperature  (°C) 

60 

Inlet  gas  composition 

Anode 

Pure  H2 

Cathode 

Air  (02:N2  =  21:79) 

Inlet  gas  flow  rate  (m3  s  1 ) 

Anode 

16.67  x  10“5 

Cathode 

25.00  x  10-5 

Thickness  of  membrane  (pan) 

30 

Size  of  catalyst  layer  (cm2) 

225 

GDL  thickness  (pan) 

300 

GDL  permeability  (m2) 

2.5  x  10“ 11 

Number  of  the  channel 

15 

Channel  width  (mm) 

1 

Shoulder  width  (mm) 

1 

Channel  depth  (mm) 

0.5,  1.0,  1.5 

Fig.  15.  Developed  view  of  single  cell  (a)  and  gas  flow  pattern  (b)  (viewpoint 
of  (b):  from  back  to  anode  separator). 


Fig.  16.  Effect  of  channel  depth  on  current  density- voltage  curve. 
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Fig.  17.  Influence  of  channel  depth  on  current  density  distribution.  Channel  depth:  (a)  0.5  mm;  (b)  1.0  mm;  (c)  1.5  mm. 


to  electrode  was  increased  by  the  gas  flow  through  GDL.  In 
this  study,  the  anode  channel  depth  was  also  changed,  how¬ 
ever,  as  the  hydrogen  concentration  overvoltage  was  ignored, 
the  gas  flow  through  anode  GDL  did  not  affect  the  reaction 
rate  directly.  And  as  the  anode  gas  of  about  100%  in  relative 
humidity  flowed  in  cell,  anode  relative  humidity  which  affect 
local  ionic  conductivity  of  membrane  became  almost  uniform 
even  if  anode  gas  flow  rate  distribution  became  remarkable 
by  the  gas  flow  through  GDL.  Consequently,  the  effect  of 
cathode  channel  depth  on  cell  performance  was  larger  than 
that  of  anode  channel  depth  under  this  operating  condition. 
Fig.  17  shows  the  current  density  distribution  of  each  chan¬ 
nel  depth  at  0.6  V.  In  this  graph,  the  current  density  was  the 
highest  at  the  turning  point  of  channel  under  all  conditions, 
and  tendency  of  0.5  mm  channel  depth  was  larger  than  other 
conditions.  However,  there  were  the  parts  of  low  current  den¬ 
sity  at  the  downstream  area.  The  reason  was  that  the  gas 
did  not  flow  uniformly  through  GDL.  It  is  considered  that 
the  non-uniform  current  density  distribution  may  cause  the 
non-uniform  water  content  and  temperature  distribution  of 
electrolyte  membrane,  and  that  it  may  reduce  the  durabil¬ 


ity  of  a  cell.  On  the  other  hand,  the  concept  of  reduction  of 
channel  depth  in  order  to  increase  the  oxygen  transfer  rate 
to  electrode  is  equal  to  that  of  the  interdigitated  channel  pro¬ 
posed  by  Nguyen  [14].  However,  the  power  of  supplying  gas 
must  be  considered  to  examine  the  effectiveness  such  sep¬ 
arator  shape.  Because  the  net  generated  power  is  different 
between  the  generated  electric  power  and  the  power  of  sup¬ 
plying.  From  these  viewpoints,  uniformity  of  current  density 
and  pressure  drop  of  each  channel  depth  were  compared. 
Fig.  18  shows  the  effect  of  channel  depth  on  the  rate  of  the 
maximum  current  density  to  the  minimum  current  density 
and  pressure  drop  at  0.6  V.  As  channel  depth  is  smaller,  the 
ratio  of  current  density  is  larger  in  Fig.  18a,  and  the  pressure 
drop  is  larger  in  Fig.  18b.  The  pressure  drop  of  the  0.5  mm 
depth  was  four  times  as  much  as  that  of  the  1.0  mm  depth. 

In  order  to  examine  the  optimal  separator  shape,  the 
performance  of  oxygen  transfer  to  electrode  and  the  pres¬ 
sure  drop  must  be  evaluated  comprehensively,  in  addition 
to  mechanical  strength,  manufacturing  cost  and  the  per¬ 
formance  of  conductivity.  Furthermore,  the  performance  of 
removing  liquid  water  in  a  channel  must  be  evaluated.  In 
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Fig.  18.  (a  and  b)  Effect  of  channel  depth  on  the  rate  of  maximum  current  density  to  minimum  current  density  and  pressure  drop  at  0.6  V. 


the  case  of  deep  channel,  it  is  expected  that  the  performance 
of  removing  liquid  water  decreases.  About  this  point,  this 
numerical  analysis  model  must  be  improved  to  examine  this. 

4.  Conclusion 

In  order  to  calculate  an  actual  scale  cell  of  PEFC  includ¬ 
ing  gas  flow  through  GDL,  the  effect  of  gas  flow  rate  through 
GDL  was  examined  by  the  numerical  analysis,  and  the  oxy¬ 
gen  transfer  model  was  made  with  the  theoretical  model. 
Next,  this  oxygen  transfer  model  was  combined  with  the 
thermal  flow  analysis  model,  and  PEFC  reaction  and  flow 
analysis  model  of  an  actual  size  was  made.  Furthermore, 
depth  of  separator  channel  was  evaluated  with  this  model 
from  the  viewpoint  of  the  characteristic  of  current  density 
and  voltage,  current  density  distribution  and  pressure  drop. 
The  validity  of  this  model  was  confirmed  by  experiments. 
The  following  results  were  obtained  by  these  examinations: 

1 .  The  oxygen  mass  transfer  rate  could  be  shown  as  the  func¬ 
tion  of  the  Reynolds  number  and  the  Schmidt  number. 
And  the  rate  was  in  promotion  of  the  square  root  of  the 
Reynolds  number. 

2.  The  experiments  without  reactions  showed  that  the  effect 
of  gas  flow  through  GDF  on  the  whole  flow  condition 
was  large  when  the  differential  pressure  between  adjoin¬ 
ing  channels  was  large. 

3.  The  numerical  analysis  showed  the  output  density 
increased  as  depth  of  separator  channel  was  smaller.  How¬ 
ever,  the  pressure  drop  and  current  density  distribution 
increased,  too. 

4.  The  current  density  at  the  turning  point  of  a  separator 
channel  was  the  highest  in  the  whole  electrode  area  under 
each  condition. 

The  numerical  analysis  developed  in  this  study  is  a  very 
simple  model  with  the  oxygen  transfer  model  and  the  one- 
dimension  flow  model  in  gas  channel.  It  does  not  suit  to 
examine  the  microscopic  phenomena,  however  the  macro¬ 
scopic  phenomena  can  be  examined,  and  this  model  can  be 


applied  to  design  of  the  PEFC  system.  On  the  other  hand, 
GDF  effective  porosity  was  used  as  a  parameter  including 
the  effect  of  liquid  water.  However,  it  is  not  advisable  to 
decide  the  GDF  effective  porosity  by  fitting  with  experimen¬ 
tal  results  under  various  conditions,  and  it  is  impossible  to 
examine  liquid  water  distribution  locally  in  GDF  from  this 
method.  Therefore,  development  of  a  detailed  liquid  water 
model  is  needed.  And  the  current  density  distribution  was 
calculated  with  the  approximate  model  of  oxygen  transfer 
in  cathode  GDF,  which  was  developed  with  the  calculation 
results  based  on  the  assumption  to  ignore  the  concentration 
diffusion  between  adjoining  channels.  (The  concentration  of 
channels  was  set  to  be  equal  to  each  other.)  For  example,  in 
the  case  of  the  turning  point  of  a  channel,  there  is  relation 
between  the  upstream  and  the  downstream  in  the  same  chan¬ 
nel,  so  the  concentration  of  the  downstream  channel  is  lower 
than  that  of  the  upstream  channel  by  electrochemical  reaction, 
therefore  it  is  necessary  to  consider  the  difference  of  concen¬ 
tration.  In  our  future  study,  it  is  expected  that  this  model 
will  be  improved  to  apply  the  various  conditions  including 
the  influence  of  liquid  water  and  various  concentration  con¬ 
ditions  in  order  to  examine  the  current  density  distribution 
strictly  in  an  actual  cell. 
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